library(tidyverse)
library(ggtext)
library(spatstat.core)
library(patchwork)
theme_set(theme_light(base_size = 10, base_family = "Open Sans"))
theme_update(
panel.grid.major = element_line(size = .3),
panel.grid.minor = element_blank(),
plot.title = element_markdown(size = 18),
plot.title.position = "plot",
plot.margin = margin(rep(25, 4))
)
Is the predictive power of the current acoustic monitoring sufficient enough to estimate the true number of killed bats?
The predictive power of the acoustic monitoring for the extrapolation of the expected number of stroke victims is decreasing
in case the crossing bats are not distributed uniformly or randomly.
with increasing length of the rotor blades.
for species with high-frequency echo calls.
We have the following parameters we can vary:
radius_rotor) - 60 m versus 33 m
area_rotor as pi * radius_rotor^2)proportion_covered)
area_monitored as area_rotor * proportion_covered)radius_monitored as sqrt(area_monitored / pi))n)
## radius circle (~ length rotor blades)
radius_rotor <- 1
## % area covered by acoustic monitoring
proportion_covered <- .05
## radius acoustic monitoring via % area covered + length rotor blades
area_rotor <- pi * radius_rotor^2
area_monitored <- area_rotor * proportion_covered
radius_monitored <- sqrt(area_monitored / pi)
## mask for window
circle <- disc(
radius = radius_rotor,
centre = c(0, 0),
mask = FALSE,
npoly = 5000
)
## number of points (~ bats)
n <- 500
from StackOverflow
sunflower <- function(n, alpha = 2, geometry = c('planar','geodesic')) {
b <- round(alpha*sqrt(n)) # number of boundary points
phi <- (sqrt(5)+1)/2 # golden ratio
r <- radius(1:n,n,b)
theta <- 1:n * ifelse(geometry[1] == 'geodesic', 360*phi, 2*pi/phi^2)
tibble(
x = r*cos(theta),
y = r*sin(theta)
)
}
radius <- function(k,n,b) {
ifelse(
k > n-b,
1,
sqrt(k-1/2) / sqrt(n-(b+1)/2)
)
}
# example:
sunflower(500, 2, 'planar')
# A tibble: 500 x 2
x y
<dbl> <dbl>
1 -0.0239 0.0219
2 0.00490 -0.0559
3 0.0440 0.0575
4 -0.0843 -0.0149
5 0.0820 -0.0521
6 -0.0279 0.104
7 -0.0538 -0.104
8 0.118 0.0430
9 -0.123 0.0509
10 0.0598 -0.128
# ... with 490 more rows
## fix random number genreration to ensure reproducibility
set.seed(12345)
## random
pp_rand <- spatstat.core::rpoint(n, win = circle)
df_rand <- as.data.frame(pp_rand) %>% mutate(dist = x^2 + y^2)
## uniform
pp_unif <- spatstat.core::runifpoint(n, win = circle)
df_unif <- as.data.frame(pp_unif) %>% mutate(dist = x^2 + y^2)
df_sunf <- sunflower(n, 1, 'planar') %>% mutate(dist = x^2 + y^2)
## directional: inner > outer
pp_inn <- rpoint(n, function(x,y) {1 - (abs(-x^2 - y^2)) + .01}, win = circle)
df_inn <- as.data.frame(pp_inn) %>% mutate(dist = x^2 + y^2)
## directional: outer > inner
pp_out <- rpoint(n, function(x,y) {(abs(x^2 + y^2)) + .01}, win = circle)
df_out <- as.data.frame(pp_out) %>% mutate(dist = x^2 + y^2)
## directional: bottom > top
pp_btm1 <- rpoint(n, function(x,y) {100 * exp(-1*y)}, win = circle)
df_btm1 <- as.data.frame(pp_btm1) %>% mutate(dist = x^2 + y^2)
pp_btm2 <- rpoint(n, function(x,y) {100 * exp(-2*y)}, win = circle)
df_btm2 <- as.data.frame(pp_btm2) %>% mutate(dist = x^2 + y^2)
pp_btm3 <- rpoint(n, function(x,y) {100 * exp(-3*y)}, win = circle)
df_btm3 <- as.data.frame(pp_btm3) %>% mutate(dist = x^2 + y^2)
pp_btm4 <- rpoint(n, function(x,y) {100 * exp(-4*y)}, win = circle)
df_btm4 <- as.data.frame(pp_btm4) %>% mutate(dist = x^2 + y^2)
pp_btm5 <- rpoint(n, function(x,y) {100 * exp(-5*y)}, win = circle)
df_btm5 <- as.data.frame(pp_btm5) %>% mutate(dist = x^2 + y^2)
pp_btm6 <- rpoint(n, function(x,y) {100 * exp(-6*y)}, win = circle)
df_btm6 <- as.data.frame(pp_btm6) %>% mutate(dist = x^2 + y^2)
pp_btm7 <- rpoint(n, function(x,y) {100 * exp(-7*y)}, win = circle)
df_btm7 <- as.data.frame(pp_btm7) %>% mutate(dist = x^2 + y^2)
pp_btm8 <- rpoint(n, function(x,y) {100 * exp(-8*y)}, win = circle)
df_btm8 <- as.data.frame(pp_btm8) %>% mutate(dist = x^2 + y^2)
pp_btm9 <- rpoint(n, function(x,y) {100 * exp(-9*y)}, win = circle)
df_btm9 <- as.data.frame(pp_btm9) %>% mutate(dist = x^2 + y^2)
pp_btm10 <- rpoint(n, function(x,y) {100 * exp(-10*y)}, win = circle)
df_btm10 <- as.data.frame(pp_btm10) %>% mutate(dist = x^2 + y^2)
plot_dist <- function(data, title, color = TRUE) {
base <-
ggplot(as.data.frame(circle), aes(x, y)) +
geom_path(color = "black", alpha = .33, size = 2) +
coord_fixed() +
scale_x_continuous(limits = c(-1, 1)) +
scale_y_continuous(limits = c(-1, 1))
if (isTRUE(color)) {
base +
geom_point(data = data, aes(color = dist), alpha = .7) +
scale_color_gradient(low = "#003F4C", high = "#B200B2", guide = "none") +
ggtitle(title)
} else {
base +
geom_point(data = data, alpha = .7) +
ggtitle(title)
}
}
## no color versions ###########################################################
## random
gg_rand_nc <- plot_dist(df_rand, "**Random**", color = FALSE)
## uniform
#gg_unif_nc <- plot_dist(df_unif, "**Uniform** (?)", color = FALSE)
gg_sunf_nc <- plot_dist(df_sunf, "**Uniform** (sunflower)", color = FALSE)
## directional: inner > outer
gg_inn_nc <- plot_dist(df_inn, "**Inner to Outer**", color = FALSE)
## directional: outer > inner
gg_out_nc <- plot_dist(df_out, "**Outer to Inner**", color = FALSE)
## directional: bottom > top
gg_btm1_nc <- plot_dist(df_btm1, "**Bottom to Top** (exp(-y))", color = FALSE)
gg_btm2_nc <- plot_dist(df_btm2, "**Bottom to Top** (exp(-2*y))", color = FALSE)
gg_btm3_nc <- plot_dist(df_btm3, "**Bottom to Top** (exp(-3*y))", color = FALSE)
gg_btm4_nc <- plot_dist(df_btm4, "**Bottom to Top** (exp(-4*y))", color = FALSE)
gg_btm5_nc <- plot_dist(df_btm5, "**Bottom to Top** (exp(-5*y))", color = FALSE)
gg_btm6_nc <- plot_dist(df_btm6, "**Bottom to Top** (exp(-6*y))", color = FALSE)
gg_btm7_nc <- plot_dist(df_btm7, "**Bottom to Top** (exp(-7*y))", color = FALSE)
gg_btm8_nc <- plot_dist(df_btm8, "**Bottom to Top** (exp(-8*y))", color = FALSE)
gg_btm9_nc <- plot_dist(df_btm9, "**Bottom to Top** (exp(-9*y))", color = FALSE)
gg_btm10_nc <- plot_dist(df_btm10, "**Bottom to Top** (exp(-10*y))", color = FALSE)
## panel
(gg_sunf_nc | gg_btm2_nc | gg_inn_nc) / (gg_rand_nc | gg_btm5_nc | gg_out_nc)
ggsave("plots/circular_distributions_nc.png", width = 15, height = 10, dpi = 700)
## colored versions ############################################################
## random
gg_rand <- plot_dist(df_rand, "**Random**")
## uniform
#gg_unif <- plot_dist(df_unif, "**Uniform** (?)")
gg_sunf <- plot_dist(df_sunf, "**Uniform** (sunflower)")
## directional: inner > outer
gg_inn <- plot_dist(df_inn, "**Inner to Outer**")
## directional: outer > inner
gg_out <- plot_dist(df_out, "**Outer to Inner**")
## directional: bottom > top
gg_btm1 <- plot_dist(df_btm1, "**Bottom to Top** (exp(-y))")
gg_btm2 <- plot_dist(df_btm2, "**Bottom to Top** (exp(-2*y))")
gg_btm3 <- plot_dist(df_btm3, "**Bottom to Top** (exp(-3*y))")
gg_btm4 <- plot_dist(df_btm4, "**Bottom to Top** (exp(-4*y))")
gg_btm5 <- plot_dist(df_btm5, "**Bottom to Top** (exp(-5*y))")
gg_btm6 <- plot_dist(df_btm6, "**Bottom to Top** (exp(-6*y))")
gg_btm7 <- plot_dist(df_btm7, "**Bottom to Top** (exp(-7*y))")
gg_btm8 <- plot_dist(df_btm8, "**Bottom to Top** (exp(-8*y))")
gg_btm9 <- plot_dist(df_btm9, "**Bottom to Top** (exp(-9*y))")
gg_btm10 <- plot_dist(df_btm10, "**Bottom to Top** (exp(-10*y))")
## panel
(gg_sunf | gg_btm2 | gg_inn) / (gg_rand | gg_btm5 | gg_out)
ggsave("plots/circular_distributions.png", width = 15, height = 10, dpi = 700)
## panel bottom > top
(gg_btm1_nc | gg_btm2_nc | gg_btm3_nc | gg_btm4_nc | gg_btm5_nc) / (gg_btm6_nc | gg_btm7_nc | gg_btm8_nc | gg_btm9_nc | gg_btm10_nc)
ggsave("plots/circular_distributions_btms_nc.png", width = 24, height = 10, dpi = 700)
## panel bottom > top
(gg_btm1 | gg_btm2 | gg_btm3 | gg_btm4 | gg_btm5) / (gg_btm6 | gg_btm7 | gg_btm8 | gg_btm9 | gg_btm10)
ggsave("plots/circular_distributions_btms.png", width = 24, height = 10, dpi = 700)
plot_monitoring <- function(gg, data, circle_monitored) {
gg +
geom_path(data = as.data.frame(circle_monitored), color = "orange2", alpha = .2, size = 1.5) +
geom_point(data = filter(data, dist < radius_monitored^2), color = "orange2")
}
plot_monitoring_all <- function(proportion_covered, df_rand, df_sunf, df_inn, df_out, df_btm_2, df_btm_5) {
## radius acoustic monitoring via % area covered + length rotor blades
area_rotor <- pi * radius_rotor^2
area_monitored <- area_rotor * proportion_covered * 2
## times 2 because we only use the lower half of the circle = half of the area later
radius_monitored <- sqrt(area_monitored / pi)
circle_monitored <- disc(
radius = radius_monitored,
centre = c(0,0),
mask = FALSE,
npoly = 5000
)
circle_half_monited <-
as.data.frame(circle_monitored) %>%
filter(y <= 0)
## random
gg_rand_r <- gg_rand_nc +
geom_path(data = circle_half_monited, color = "orange2", alpha = .2, size = 1.5) +
geom_point(data = filter(df_rand, dist < radius_monitored^2 & y <= 0), color = "orange2")
## uniform
# gg_unif_r <- gg_unif_nc +
# geom_path(data = as.data.frame(circle_monitored), color = "orange2", alpha = .2, size = 1.5) +
# geom_point(data = filter(df_unif, dist < radius_monitored^2), color = "orange2")
gg_sunf_r <- gg_sunf_nc +
geom_path(data = as.data.frame(circle_monitored), color = "orange2", alpha = .2, size = 1.5) +
geom_point(data = filter(df_sunf, dist < radius_monitored^2), color = "orange2")
## directional: inner > outer
gg_inn_r <- gg_inn_nc +
geom_path(data = as.data.frame(circle_monitored), color = "orange2", alpha = .2, size = 1.5) +
geom_point(data = filter(df_inn, dist < radius_monitored^2), color = "orange2")
## directional: outer > inner
gg_out_r <- gg_out_nc +
geom_path(data = as.data.frame(circle_monitored), color = "orange2", alpha = .2, size = 1.5) +
geom_point(data = filter(df_out, dist < radius_monitored^2), color = "orange2")
## bottom > top
gg_btm2_r <- gg_btm2_nc +
geom_path(data = as.data.frame(circle_monitored), color = "orange2", alpha = .2, size = 1.5) +
geom_point(data = filter(df_btm2, dist < radius_monitored^2), color = "orange2")
gg_btm5_r <- gg_btm5_nc +
geom_path(data = as.data.frame(circle_monitored), color = "orange2", alpha = .2, size = 1.5) +
geom_point(data = filter(df_btm5, dist < radius_monitored^2), color = "orange2")
## panel
p <- (gg_sunf_r | gg_btm2_r | gg_inn_r) / (gg_rand_r | gg_btm5_r | gg_out_r)
}
(panel_monitored <- plot_monitoring_all(
proportion_covered = .5, df_rand, df_sunf, df_inn, df_out, df_btm_2, df_btm_5
))
ggsave("plots/circular_distributions_area_monitored_50.png", width = 15, height = 10, dpi = 700)
(panel_monitored <- plot_monitoring_all(
proportion_covered = .4, df_rand, df_sunf, df_inn, df_out, df_btm_2, df_btm_5
))
ggsave("plots/circular_distributions_area_monitored_40.png", width = 15, height = 10, dpi = 700)
(panel_monitored <- plot_monitoring_all(
proportion_covered = .3, df_rand, df_sunf, df_inn, df_out, df_btm_2, df_btm_5
))
ggsave("plots/circular_distributions_area_monitored_30.png", width = 15, height = 10, dpi = 700)
(panel_monitored <- plot_monitoring_all(
proportion_covered = .2, df_rand, df_sunf, df_inn, df_out, df_btm_2, df_btm_5
))
ggsave("plots/circular_distributions_area_monitored_20.png", width = 15, height = 10, dpi = 700)
(panel_monitored <- plot_monitoring_all(
proportion_covered = .1, df_rand, df_sunf, df_inn, df_out, df_btm_2, df_btm_5
))
ggsave("plots/circular_distributions_area_monitored_10.png", width = 15, height = 10, dpi = 700)
(panel_monitored <- plot_monitoring_all(
proportion_covered = .05, df_rand, df_sunf, df_inn, df_out, df_btm_2, df_btm_5
))
ggsave("plots/circular_distributions_area_monitored_05.png", width = 15, height = 10, dpi = 700)
Sys.time()
[1] "2021-11-22 14:31:49 CET"
R version 4.1.0 (2021-05-18)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 10 x64 (build 19043)
Matrix products: default
locale:
[1] LC_COLLATE=German_Germany.1252 LC_CTYPE=German_Germany.1252
[3] LC_MONETARY=German_Germany.1252 LC_NUMERIC=C
[5] LC_TIME=German_Germany.1252
system code page: 65001
attached base packages:
[1] stats graphics grDevices utils datasets methods
[7] base
other attached packages:
[1] patchwork_1.1.1 spatstat.core_2.3-0 rpart_4.1-15
[4] nlme_3.1-152 spatstat.geom_2.3-0 spatstat.data_2.1-0
[7] ggtext_0.1.1 forcats_0.5.1 stringr_1.4.0
[10] dplyr_1.0.7 purrr_0.3.4 readr_2.0.2
[13] tidyr_1.1.4 tibble_3.1.5 ggplot2_3.3.5
[16] tidyverse_1.3.1
loaded via a namespace (and not attached):
[1] fs_1.5.0 spatstat.sparse_2.0-0 lubridate_1.8.0
[4] httr_1.4.2 tools_4.1.0 backports_1.2.1
[7] bslib_0.3.1 utf8_1.2.2 R6_2.5.1
[10] DBI_1.1.1 mgcv_1.8-35 colorspace_2.0-2
[13] withr_2.4.2 tidyselect_1.1.1 downlit_0.2.1
[16] compiler_4.1.0 textshaping_0.3.6 cli_3.0.0
[19] rvest_1.0.2 xml2_1.3.2 labeling_0.4.2
[22] sass_0.4.0 scales_1.1.1 goftest_1.2-3
[25] systemfonts_1.0.3 digest_0.6.28 spatstat.utils_2.2-0
[28] rmarkdown_2.11 pkgconfig_2.0.3 htmltools_0.5.2
[31] highr_0.9 dbplyr_2.1.1 fastmap_1.1.0
[34] rlang_0.4.12 readxl_1.3.1 rstudioapi_0.13
[37] farver_2.1.0 jquerylib_0.1.4 generics_0.1.0
[40] jsonlite_1.7.2 distill_1.3 magrittr_2.0.1
[43] Matrix_1.3-3 Rcpp_1.0.7 munsell_0.5.0
[46] fansi_0.5.0 abind_1.4-5 lifecycle_1.0.1
[49] stringi_1.7.5 yaml_2.2.1 grid_4.1.0
[52] crayon_1.4.1 deldir_1.0-5 lattice_0.20-44
[55] haven_2.4.3 splines_4.1.0 tensor_1.5
[58] gridtext_0.1.4 hms_1.1.1 knitr_1.36
[61] pillar_1.6.4 markdown_1.1 reprex_2.0.1
[64] glue_1.4.2 evaluate_0.14 modelr_0.1.8
[67] vctrs_0.3.8 tzdb_0.1.2 cellranger_1.1.0
[70] gtable_0.3.0 polyclip_1.10-0 assertthat_0.2.1
[73] xfun_0.27 broom_0.7.9 ragg_1.1.3
[76] ellipsis_0.3.2